分子系统发育学对传染病的研究产生了深远的影响,特别是对RNA病毒等快速进化的传染病病原体的研究。它洞察了流行病暴发和季节性疾病的起源、进化史、传播途径和来源人群。关于快速进化病毒的一个关键观察是,进化和生态过程发生在相同的时间尺度上(Pybus和Rambaut, 2009)。这一点之所以重要,有两个原因。首先,它意味着中性遗传变异可以跟踪生态过程和种群动态,提供过去的进化事件(如系谱关系)和过去的生态/种群事件(地理传播和种群规模和结构的变化)的记录,这些事件没有直接观察到。第二,进化和生态过程的同时导致了它们的相互作用,当它们不是平凡的时候,就需要联合分析。
可以说,迄今为止研究最多的传染病病原体是人体免疫缺陷病毒(艾滋病毒),它是数以千计的系统发育研究的主题。这些研究揭示了艾滋病毒进化生物学、流行病学、起源、系统地理学、传播动力学和耐药性的许多方面。事实上,大量关于艾滋病毒的文献表明,在病毒进化的背景下,可以更好地理解一种迅速进化的病原体的生物学的几乎每一个方面。无论是追溯艾滋病毒大流行的人畜共患病起源,还是描述病毒种群和宿主免疫系统之间的相互作用,系统发育分析常常能提供线索。
尽管系统发育的概率建模方法早于Sanger测序(Edwards and Cavalli-Sforza, 1965),但直到最近十年,概率建模才成为系统发育重建的主要方法。这种主导地位的部分原因是贝叶斯推理(Huelsenbeck等人,2001年)的兴起,它在描述先验知识方面具有很大的灵活性,能够通过Metropolis-Hastings算法应用于复杂的高参数模型,并且可以轻松地将多个数据源集成到单一分析中。分子进化和系统发育的概率模型的历史是一个逐渐完善的历史;选择那些在描述不断增长的经验数据方面具有最大效用的建模变量的过程。一个新模型的效用要么通过它与数据的拟合程度(正式模型比较或拟合优度测试)进行评估,要么通过它允许研究人员对数据提出的新问题进行评估。在这篇综述中,我们将描述传染病领域的现代系统发育方法,特别是参考快速进化的病毒病原体如丙型肝炎病毒(HCV)、艾滋病毒和甲型流感病毒的系统发育流行病学的贝叶斯推断。本综述分为两个主要部分。在第2节中,我们讨论了重建传染病历史的系统发育方法,包括起源的确定,共同祖先的年代,宽松的系统发育和基于合并的种群动态。在第3节中,我们回顾了流行病学模型,最后概述了将统计系统发育与动力学建模相结合的系统动力学模型的发展进展。
引入一种有效的方法,在给定的系统发育树中计算序列对齐的概率(称为系统发育似然;Felsenstein, 1981)预示着统计框架下实际系统发生树重建的开始。大约与此同时,“聚合”理论也被引入:一种将随机样本个体的族谱形状与他们所来自的人口规模联系起来的理论(Kingman, 1982;详情见第2.3节)。这两项进展后来发展到这样的地步,它们一起使估计病毒的进化历史和过去的种群动态成为可能。
贝叶斯推理将似然Pr(D∣θ)(给定模型参数时数据的概率)和先验P(θ)(在看到数据之前模型参数的概率)集合在一起,因此给定数据时模型参数的后验概率(θ)为: $$ P(\theta \mid D)=\frac{\operatorname{Pr}(D \mid \theta) P(\theta)}{\int \operatorname{Pr}(D \mid \theta) P(\theta) d \theta} $$
在标准的系统发育设置中,概率模型参数包括系统发育树、聚合时间和替代参数,必须指定这些参数的先验概率分布。通过使用Kingman 's coalescent作为树的先验密度,贝叶斯推断可以用来同时估计病毒序列的系统发育和病毒种群的人口统计历史(Drummond et al., 2002, Drummond et al., 2005, Drummond et al., 2006,见框1)。系统发育推断方法的扩展以适应时间标记的序列数据(Rambaut, 2000, Drummond et al., 2002)和放宽严格分子钟的假设(Thorne et al., 1998,Kishino et al., 2001, Sanderson, 2002, Drummond et al., 2006, Rannala and Yang, 2007)提供了复杂的祖先散度时间估计方法。对于占据多个宿主物种的病毒种(如甲型流感),旨在检测跨物种传播的模型可能为宿主种群中病毒株的起源提供线索(Reis等人,2009年)。
当一种新的流行病出现时,首要目标之一是追溯其基因和地理起源。重建系统发育树以推断进化关系一直是发现区域流行病起源的关键工具,例如由HIV (Gao等人,1999,Santiago等人,2002)、HCV (Pybus等人,2009,Markov等人,2009)和SARS冠状病毒(SARS- cov) (Li等人,2005)引起的流行病。一些研究也试图使用系统发育树来得出关于病毒流行的传播历史和地理传播的结论(Motomura等人,2003年,Santiago等人,2005年,Gilbert等人,2007年)。然而,在对未在系统发育树重建中明确建模的流行病过程方面作出结论时,应非常小心,即使是这样,用户也需要考虑基本模型假设的适当性。
用于确定流行病起源的一种常见而直接的方法是确定与流行病最密切相关的非流行病基因型或谱系,即在系统发育树中与流行病毒株聚集最紧密的分子序列。虽然这种方法是直观的,但它的成功很大程度上取决于收集到的数据。
与HIV-1最接近的猿类免疫缺陷病毒(SIV)是SIVcpz (Gao等人,1999,Santiago等人,2002),它存在于黑猩猩亚种Pan troglodytes troglodytes和P.t. schweinfurthii中,以各自亚种特异性SIV谱系SIVcpzPtt和SIVcpzPts的形式存在。尽管SIVcpz在被确认后就成为HIV-1人畜共患病来源的主要候选,但由于黑猩猩感染的数量较少,不能排除其他来源(Vanden Haesevelde等人,1996年)。很久以后,从喀麦隆森林中野生P. t. troglodytes猿类的粪便样本中收集到SIVcpz后,才确认了HIV-1的来源(Keele et al., 2006)。HIV-1组M和N与来自粪便样本的序列的关系比之前确定的SIVcpz菌株密切得多。这一发现揭示了HIV-1 M组(大流行)和N组(非大流行)的不同起源,分别溯源于喀麦隆东南部和中部的黑猩猩群落。通过系统发育技术对这些野生黑猩猩HIV-1宿主的精确地理识别,为SIVcpz导致艾滋病毒/艾滋病流行提供了关键证据。
相反,如果不能确定与流行毒株足够密切相关的毒株,那么系统发育树就不能轻易地提供关于起源的答案。例如,关于1918年甲型H1N1流感大流行的起源——它的来源是禽类、非人类哺乳动物,甚至是人类——有很多激烈的争论。这种不确定性主要源于1918年病毒的直接祖先源种群缺乏序列(Gibbs and Gibbs, 2006)。
在寻找HIV-1 O组的起源时也遇到了类似的问题,尽管没有那么严重。研究发现,HIV-1 O组毒株与西部低地大猩猩(大猩猩大猩猩)中的SIVgor关系最为密切(Van Heuverswyn et al., 2006, Takehisa et al., 2009)。然而,HIV-O序列与已知的SIVgor序列有一定的差异,因此,导致HIV-1 O组和SIVgor的传播途径仍不确定。
在系统发育树中,一个新出现的病毒株与其他病毒株的穿插通常被解释为支持多个独立病毒传入的证据。例如,HIV谱系是副系的,而SIV谱系产生了几个独立的HIV集群,表明多重人畜共患病毒传播到人群中(Santiago等人,2005年,Keele等人,2006年)。虽然从直觉上看,突发病毒的单独群集表明有多次传入,但仅从群集的数量还不清楚有多少独立事件导致了观察到的模式。不完整的分类单元采样将导致少计数。例如,可能存在一个未采样的序列将分裂一个紧急病毒集群,或一个额外的未采样的紧急病毒集群。这两种情况,如果检测到,将增加推断事件数量的下界。由于系统发育估计误差,事件的数量也可能被错误估计。最后,在事件可能可逆的情况下,如耐药突变,例如H3N2流感病毒中的金刚烷耐药(Nelson et al., 2009),很可能在系统发育史中也存在逆转,而这些逆转并不总是通过简单的简约性重建检测到,再次导致少计数。基于所有这些原因,在系统发育上应用贝叶斯系统地理建模和特征进化对定量评估这些不同误差来源产生的不确定性至关重要(见第2.4节)。
与HIV-1相比,近二十年来已经明确确定,HIV-2的前身是来自黑白眉(Cerocebus torquatus atys)的SIVsm (Hirsch et al., 1989, Gao et al., 1992)。(Santiago et al., 2005)根据HIV-2株的系统发育树和分支位置显示的清晰的地理聚类,认为HIV-2 A组和B组的地理起源在东部黑白眉范围。尽管这种启发式方法被普遍用于定位系统地理起源,但除了前面提到的抽样误差外,它还有几个缺点。首先,它依靠强大的地理信号在树木中产生明确的地理聚集模式。其次,缺乏正式的统计框架导致无法量化与地理估计相关的不确定性。许多统计系统发育方法旨在通过将地理位置作为沿着树进化的另一个状态来重建迁移过程。状态要么是离散的(Lemey et al., 2009b),用城市或省份的名称表示,要么是连续的,用地点的经纬度表示(Biek et al., 2006, Biek et al., 2007, Lemey et al., 2010)。
即使是全面采样,单一的系统发育树也不足以反映经过重组或重组的病毒种的复杂遗传起源。当病毒基因组的片段来自不同的病毒时,就会发生重组,而重组也需要来自一个来源的遗传物质与来自另一个来源的遗传物质(断开)结合。这两个过程使从两种现有基因型产生新的组合成为可能。此外,这些通常较大的基因变化可能为适应新的宿主物种提供了潜力(Parrish et al., 2008)。重组在A型流感病毒的进化中发挥了重要作用(Lindstrom等人,2004年,Holmes等人,2005年,Nelson等人,2008年)。在登革热(Holmes等人,1999年)、艾滋病毒、HCV和SARS-CoV中也发现了重组的证据(Li等人,2006年)。
有许多系统发育方法旨在通过识别排列不同部分的拓扑不一致来检测重组(Grassly和Holmes, 1997, Salminen等人,1995,Lole等人,1999,Smith, 1992, Robertson等人,1995,Paraskevis等人,2005),这是重组的一个潜在后果。这些方法大多使用滑动窗口方法沿序列长度计算汇总统计信息。系统发育方法是基于估计(i)每个窗口的自举值或(ii)枝后验概率,自举值的突然变化,枝后验概率或位点百分比标识是区域周围存在断点的指示。其他方法显式估计断点在对齐中的位置,提供了测试对重组支持强度的访问(Holmes等人,1999年)。最后,一些方法通过网络来描述进化历史,包括水平转移(Huson, 1998)或祖先重组图(Bloomquist和Suchard, 2010)。
一般来说,RNA病毒突变迅速,因此仅相隔几个月分离的病毒可能表现出可测量的遗传差异(Drummond等人,2003a和其中的参考文献)。事实上,一些RNA病毒的突变率是如此之高,以至于在感染过程中会导致宿主内部的进化变化。艾滋病毒和丙型肝炎病毒等病毒引起的长期慢性感染尤其如此。因此,对相隔数年采样的序列进行分析,就好像它们是同时期的一样,是不适当的。具有这种时间结构的序列数据被称为异时数据,从这样的数据可以估计替代率和发散时间校准到日历尺度。在这里,以日历单位表示分支长度的树称为“时间树”。图1描述了快速进化病毒的串行采样时间树示例。
为了考虑序列数据中的时间结构,最早的方法是通过估计一个基因树的分枝长度不受限制,然后对根到尖端的遗传距离与采样时间进行线性回归来估计时间尺度(见Drummond等人,2003b)。该方法用于首次估计HIV-1 M组的最近共同祖先(tMRCA)发生的时间,将其确定在20世纪30年代(Korber et al., 2000)。尽管这种方法很简单,但它也准确地估计了1959年样本中最古老的HIV序列的年龄。一种基于极大似然的方法(单速率日期tips (SRDT)模型;Rambaut, 2000),假设一个严格的分子钟,估计了一个固定树上的祖先分化时间和总体替代率。SRDT模型用于确定HIV-2亚型A的最近共同祖先为1940±16年,B亚型为1945±14年(Lemey et al., 2003年)。在贝叶斯合并方法(Drummond et al., 2002, Drummond et al., 2005, Drummond and Rambaut, 2007)中使用序列合并作为树先验,允许时间尺度与其他系统发育和人口统计参数同时估计。最近,一种宽松的时钟贝叶斯凝聚分析,包括1959年(ZR59)和1960年(DRC60)的两个历史病毒样本(Worobey等人,2008年),将HIV-1 M组的tMRCA估计推迟到1908年(1884-1924)。
除了估计疫情暴发的时间外,了解疫情毒株的祖先在疫情暴发前在源人群中传播了多长时间也可能很重要。这有时可以通过流行支的祖先支的长度来表示。以2009年猪源型A型流感病毒为例,根据所分析的病毒片段,导致2009年S-OIV毒株的分支长度估计为9-17年,这表明未采样的多样性大约为10年(Smith et al., 2009)。
为了估计SIVsm株共同祖先的年龄,对HIV-2/SIVsm的tMRCA进行了年代测定,表明在HIV-2 A组和B组人畜共患病之前的共同祖先仅跨越了过去几个世纪(Wertheim和Worobey, 2009)。这并不一定表明SIVsm最早出现在几个世纪以前,只是表明所有当前SIVsm的共同祖先可能是最近才出现的。然而,就连这一结论最近也受到了质疑(Worobey等人,2010年),因为独立校准证据表明tMRCA实际上可能大于32000年前,这导致了对通常用于散度时间定年的统计替代模型的保真度的争论,因为真实的散度时间与采样间隔相比非常古老。Wertheim和Pond在2011年证明,不考虑选择的影响的替代模型会产生低估的分支长度,导致在存在净化选择的情况下,年龄估计要小得多。如果数据集的总采样间隔只是树的总年龄的一小部分,那么问题就更大了。
虽然纳入采样日期为系统发育推断提供了额外的信息,但这也意味着这些日期的可靠性对推断的有效性有很大的影响。1977年重新出现的H1N1流感病毒被发现错过了几十年的进化,在基因上与1950年的H1N1病毒非常相似(Nakajima et al., 1978)。因此,它被认为是一种在未知实验室中冷冻了几十年之后再次成为“野生”菌株的后裔(Zimmer和Burke, 2009)。如果没有对缺失的进化进行修正,包括再出现菌株在内的分析将产生有偏差的日期估计,并增加再出现谱系和整个系统发育的tMRCA的方差(Wertheim, 2010)。在序列的采样日期有争议或未知的情况下,需要一种可以处理日期未知的序列的方法。例如,叶年代测定法将序列的未知日期或年龄作为一个参数,将其与内部节点的年龄处理相同(Drummond et al., 2003c, Nicholls and Gray, 2008, Shapiro et al., 2010)。
不现实的采样日期也可能是人为错误的结果,因此在分析之前没有被识别。因此,对不现实的日期进行诊断对于发现记录日期中的错误非常重要。一种可能的方法是,如果病毒没有显示出显著偏离恒定速率(Wertheim, 2010),根据采样年绘制根到尖的遗传距离图。另一种方法是通过依次降低每个校准点并重新估计日期来检查校准,以确认估计的日期是一致的(Shapiro et al., 2010, Ryder and Nicholls, 2011)。
容纳异时数据的早期方法假定了一个严格的时钟模型。然而,一项使用SRDT模型对异时RNA病毒序列的综合研究(Rambaut, 2000)表明,在研究的50种RNA病毒中,大多数都拒绝恒定速率分子时钟假说(Jenkins et al., 2002)。无根系统发育是系统发育树各分支的速率变异性的另一个极端。它们都不是潜在进化过程的现实表现,而现实则介于两者之间。这催生了许多方法的发展,这些方法放松了分子钟的假设,并在它们对不同分支的速率变化模式的假设上有所不同。
本地时钟模型方法为树的分支/区域分配不同的速率。然而,在没有外部信息的情况下,很难先验地知道将树划分为本地时钟模型的最佳方法。贝叶斯模型平均通过对所有可能的本地时钟模型(Drummond和Suchard, 2010)进行平均,克服了速率分配的挑战,同时估计替代率,以及替代率变化的数量和位置。
另一类放松时钟模型基于“速率平滑”,包括非参数速率平滑(Sanderson, 1997)、惩罚似然(Sanderson, 2002)和贝叶斯自相关放松时钟方法(Thorne等人,1998,Kishino等人,2001,Aris-Brosou和Yang, 2002, Rannala和Yang, 2007)。这些方法通过惩罚与父分支速率较大的偏差来限制父分支和子分支的速率相似。因此,速率的变化是通过小而频繁的变化发生的。不同的贝叶斯自相关时钟模型在给定父速率的情况下,用于建模分支速率的分布是不同的(Thorne等人,1998年,Kishino等人,2001年)。
然而,对甲型流感和登革热4型的序列数据的分析并没有提供分支率自相关的任何证据(Drummond等人,2006年),这表明自相关模型可能不适用于分析来自单一病毒种的序列谱系。尽管谱系效应可能会导致速率的自相关(通过生活史的渐进变化,代谢速率等),但达尔文选择的基因特异性作用也会在谱系之间引起明显的速率变化,通过在整个系统发育过程中产生分子钟的普遍过度分散(Takahata, 1987, Takahata, 1991)。谱系之间速率变化的第二个来源可以用不相关的放松时钟模型更好地模拟(Drummond等人,2006年),该模型没有对祖先和后代分支之间速率的自相关做假设。已发表的分析提供了强有力的证据,支持不相关的放松时钟模型(例如,Salemi等人,2008年,Worobey等人,2008年)胜过严格时钟模型。
除了估计祖先差异的年龄外,如果疾病起源于人畜共患病,估计跨物种传播的时间也很有意义。确定主机切换时间的一种方法是应用非齐次替换模型。非同质替代模型的动机是承认病毒在不同宿主物种内的替代模式可能存在差异,这违反了标准替代模型基础上的同质性和平稳性假设。因此,对树的不同部分应用不同的替代模型可能更合适(Forsberg和Christiansen, 2003)。非齐次替代模型允许平衡频率,因此模型参数在一个分支上发生变化,并且从变化点开始的所有后代谱系都假设在该点之前的谱系具有不同的平衡基频率。该技术已被用于表明1918年A型流感病毒的直系祖先种群居住在哺乳动物宿主中(Reis等人,2009年)。然而,它并没有表明猪流感病毒和1918年流感病毒的最近共同祖先是否存在于人类或其他哺乳动物体内。
解释估计的发散时间可能很困难。可能有更古老的直系祖先,但由于基因漂变等过程,可以揭示他们的血统没有被采样或没有存活到现在。因此,估计的tMRCA可能不能回答兴趣的问题。对于由人畜共患病传播引起的流行病,宿主切换事件是最重要的,但估计流行毒株的tMRCA不能直接估计传播时间,而只是作为一个下界。同样,如果过去曾有过导致遗传多样性丧失的过程,或者采样不全面,那么估计的tMRCA可能大大小于病毒谱系的年龄。前者的一个明显例子发生在季节性流感中,这是由于季节性人口波动造成的,也包括免疫监视引起的强正达尔文选择(Fitch等人,1991年,Bush等人,1999年),导致任何单季节样本的快速谱系更替和最近的共同祖先。
类似地,Worobey等人(2008)的分析表明,HIV-1 M组的tMRCA似乎被推迟了,因为加入了一个来自1960年的额外的流行前样本,该样本与1959年的序列(ZR59)高度不同。一般来说,包含较老的样本可以增加估计的根的年龄,因为(i)揭示了以前未采样的谱系,这些谱系是没有它们的tMRCA估计的外群,或者(ii)只是因为更多的时间采样破坏了较长的内部分支,并可能揭示了被认为是现代的变体的古老证据,导致估计速率较慢,因此估计的根高较老。
最后,单靠目前的技术很可能不能总是恢复遥远过去的准确的分歧日期,正如最近的分析所表明的那样,SIV的历史比以前认为的要深刻得多(Worobey等人,2010)(Sharp等人,2000,Wertheim和Worobey, 2009)。图2用三棵估计的病毒时间树说明了这个问题,它们的最近共同祖先的推断年龄相差很大。我们期望对甲型流感时间树的推断年龄有最大的置信度,因为样本周期是时间树总年龄的很大一部分,而对丙型肝炎时间树的推断年龄的置信度最低,因为样本周期是大于1000年的推断年龄的一小部分。
因此,除了更好的跨谱系变异性速率模型(参见Guindon等人,2004年,这一方向的早期步骤),未来对分化时间定年的研究可能将关注于更准确地解释纯化选择及其在维持编码基因的结构和功能方面的作用的模型。达尔文选择的影响既表现为谱系的扭曲(O’fallon, 2010, O’fallon等人,2010),也表现为中性预期的替代过程(例如,Bloom等人,2007,Cartwright等人,2011)。在易受克隆干扰、结构紧凑、信息丰富、进化轨迹受功能和结构高度约束的病毒基因组中,特别是考虑长时间周期时,考虑普适纯化选择的作用尤为重要。除此之外,还需要采用统计上更为严格的方法,以纳入各种校准信息来源,例如生物地理学、考古学和古生物学证据。贝叶斯统计框架特别适合这种多信息源的集成。
基于家谱的种群遗传学可以用来推断人口参数,包括种群大小、增长或下降的速度和种群结构。当人口波动的特征时间尺度与取代的累积率相当时,过去的人口动态就被“记录”在分子序列的取代模式中。因此,聚合理论可以与异时序列的时间信息相结合,以发现过去的流行病学事件,并在日历时间尺度上确定它们。
Kingman 's coalescent (Kingman, 1982)描述了假设理想化的Wright - Fisher种群(Fisher, 1930, Wright, 1931),样本谱系中的coalescent时间与种群规模之间的关系。最初的公式适用于恒定的总体,但该理论已推广到任何确定性变化的总体大小函数,其积分$\int_{t_0}^{t_1} N(t)^{-1}$计算(格里菲斯和Tavaré, 1994)。具有预先定义的种群函数的参数化模型,如指数增长、扩展模型和逻辑增长模型,可以很容易地用于一个合并框架中(详见图3和框1)。例如,在贝叶斯凝聚框架中使用“分段逻辑”人群模型来估计埃及HCV基因4a型感染的人群历史(Pybus等,2003年)。这一分析表明,在1930-1955年期间,丙型肝炎病毒在埃及迅速蔓延,这与实施抗血吸虫病注射的公共卫生运动导致丙型肝炎病毒在埃及蔓延的假设是一致的。
聚合过程是高度可变的,因此对多个不相关的位点进行采样(Felsenstein, 2006, Heled和Drummond, 2008)或增加采样时间的时间分布(Seo等人,2002)都可以用来提高基于聚合的方法的统计能力,并提高估计种群大小和替代率的精度(Seo等人,2002)。然而,在许多病毒种类中,整个基因组作为一个位点,或只有在通过重复感染出现机会时才进行重组。因此,由于缺乏独立的位点,对种群历史估计的精度有了一个上限。
在许多情况下,种群规模历史的精确函数形式是未知的,简单的种群增长函数可能不能充分描述相关的种群历史。非参数聚合法通过直接从序列数据估计作为时间函数的种群大小提供了更大的灵活性,并且可以用于数据探索,以指导进一步分析的参数种群模型的选择。这些方法首先将时间树分割成段,然后根据每个段之间的合并区间分别估计每个段的种群大小。
这些方法之间的主要区别是(i)种群规模函数如何沿树分割,(ii)所使用的统计估计技术,以及(iii)在贝叶斯方法中,控制种群规模函数的参数的先验密度的形式。在“经典天际线图”(Pybus et al., 2000)中,每个合并区间被视为一个单独的段,因此n个类群的树具有n−1个种群大小参数。然而,种群大小变化的真实数量可能要少得多,广义天际线图(strmmer和Pybus, 2001)通过根据小样本Akaike信息准则(AICc)对区间进行分组(Burnham和Anderson, 2002)承认了这一点。使用广义天际线图(Strimmer和Pybus, 2001)调查了HIV-2的流行史,表明在几内亚比绍HIV-2亚型A的早期历史中,人口规模相对稳定,直到最近才扩大(Lemey等人,2003)。利用这些信息,作者然后采用分段扩展增长模型,估计扩展到1955-1970年范围的时间。
虽然广义天际线图是数据探索的好工具,并有助于模型选择(例如,Pybus等人,2003年,Lemey等人,2004年),但它基于单一输入树推断人口历史,因此不考虑系统发育重建产生的抽样误差,也不考虑合并过程的内在随机性。通过在贝叶斯统计框架中实现天际线图方法,同时推断样本谱系、替代参数和种群规模历史,克服了这一缺点。广义天际线图的进一步扩展包括用分段线性函数而不是分段常数总体来建模总体规模,允许随着时间的持续变化而不是突然的跳跃。贝叶斯天际线图(Drummond等人,2005年)曾被用来表明,HIV-1 M组的有效人群规模可能在20世纪上半叶以相对较慢的速度增长,随后是更快的增长(Worobey等人,2008年)。在更短的时间范围内,对从一对HIV-1供体和受体收集的数据集进行贝叶斯天际线图分析,揭示了病毒传播后基因多样性的大量丧失(Edwards等人,2006年)。用恒定逻辑增长模型进一步分析估计,在水平传播过程中,供体中存在的HIV-1基因多样性损失了99%以上。这具有重要的意义,因为瓶颈下的过程决定了病毒在受体宿主中的适应性。
贝叶斯天际线图的一个缺点是,人口规模的变化数量必须由用户先验地指定,而适当的数字很少为人所知。一种解决方案是利用可逆跳跃MCMC (Opgen-Rhein等人,2005年)或贝叶斯变量选择(Heled和Drummond, 2008年)对人口模型进行贝叶斯模型平均,在这种情况下,人口规模变化的数量是作为模型一部分估计的随机变量。
到目前为止讨论的人口统计推断方法都假定在感兴趣的人口中没有细分。与大小的变化一样,种群结构也会对聚结区间大小的模式产生影响,因此,当种群结构存在时,结果的可靠性会受到质疑(Pybus et al., 2009)。
系统地理学是研究种群或类群空间分布的演化和扩散过程的学科。系统地理学方法可分为两种途径。第一种方法进行树重建后分析以回答系统地理问题,第二种方法联合估计感兴趣的系统发育和系统地理参数。当将地理位置作为离散状态处理时,前一种方法在过去几十年里很受欢迎。它的优点是计算量较少,但分析的结果取决于输入树。由于其简单性,推断祖先位置最流行的方法是最大限度简约(Slatkin和Maddison, 1989, Swofford, 2003, Maddison和Maddison, 2005, Wallace等人,2007),然而,这种方法不允许对与祖先位置重建相关的不确定性进行任何概率评估。
迁移模型是用于分析迁移过程的突变模型。最近一项关于甲型H5N1流感病毒的研究引入了一种全概率“迁移”方法,通过连续时间马尔可夫过程对病毒谱系的地理移动过程进行建模,其中状态空间由序列采样的位置组成(Lemey等人,2009b)。这有助于估计成对位置之间的迁移率。此外,该方法估计了树中内部节点的祖先位置,并使用贝叶斯变量选择(BVS)推断主要迁移路线,并对不同位置(或宿主群体)之间的连通性不确定性提供模型平均。这种方法有助于调查甲型H5N1流感的起源及其全球传播途径,也有助于重建新型甲型H1N1流感大流行的最初传播(Lemey et al., 2009b)。然而,离散位置状态模型的一个共同局限性是祖先位置仅限于采样位置。对西非和中非犬狂犬病数据集的分析表明,缺乏靠近根部的序列采样可能会阻碍对病毒地理起源的准确估计(Lemey等人,2009b)。因此,通过增加采样的空间密度和时间深度,系统地理估计得到了改进。然而,密集的地理采样导致了大的系统发育和计算密集的分析。
结构化聚合(Hudson, 1990)也可用于研究系统地理学。结构化聚合也被扩展到异时数据(Ewing等人,2004年),从而可以估计日历单位中demes之间的迁移率。序列结构化聚合第一次应用于具有两个底体的HIV数据集,以研究患者内亚群体的动态(Ewing等人,2004年),但同样类型的推断也可以在宿主群体水平上进行。模型的进一步发展允许demes的数量随时间变化(Ewing和Rodrigo, 2006a)。MIGRATE (Beerli和Felsenstein, 2001)还使用结构化聚并法在贝叶斯和极大似然框架下估计亚种群规模和迁移率,最近被用于调查病毒流行的空间特征(Bedford等人,2010)。此外,一些研究聚焦于幽灵demes的影响(Beerli, 2004, Ewing和Rodrigo, 2006b),但尚无明确包含种群结构、异时样本和非参数种群规模历史的模型。
一种特别的解决方案涉及沿着树的迁移过程建模,该方法有条件地独立于由天际线图估计的种群规模(Lemey et al., 2009a)。因此,给定树,迁移过程被认为与合并先验无关。然而,这种方法并没有捕捉到结构性融合中隐含的迁移和融合之间的相互作用,因为融合率应该取决于谱系所处的demme的种群规模。正如我们将在下一节中看到的,统计系统地理学是一个将系统发育和数学流行病学模型统一起来看起来非常有前途的领域。
在某些情况下,将样本的空间方面建模为连续变量更为合适。野生动物宿主种群的系统地理学通常是通过扩散模型在空间连续体中建模的,因为病毒传播和宿主移动往往难以用少量离散的demes来模拟。一个例子是浣熊特异性狂犬病毒在美国东部的地理范围扩大(Biek等人,2007年,Lemey等人,2010年)。通过比较方法(Felsenstein, 1985年,Harvey和Pagel, 1991年),布朗扩散也被用于模拟从蒙大拿州西部美洲狮(Puma concolor)种群中收集的猫免疫缺陷病毒的系统地理。由于病毒主要垂直传播,由此产生的系统地理重建被用作宿主人口历史和种群结构的代理(Biek等人,2006年)。然而,布朗扩散的一个假设是所有分支上的速率同质性。通过将放松时钟模型的概念扩展到扩散过程,可以放宽这一假设(Lemey等人,2010年)。仿真结果表明,当空间运动的基本过程类似于过分散的随机游走时,松弛扩散模型比布朗扩散模型具有更好的覆盖率和统计效率。
就像他们的迁移模型一样,这些模型在形成样本谱系时忽略了人口密度和地理分布的相互作用。然而,将合并框架扩展到空间连续体的数学理论的发展已经取得了进展(Barton等人,2002,Barton等人,2010a, Barton等人,2010b),尽管还没有开发出在这些模型下提供推理的方法。
用MCMC分析贝叶斯聚结分析. 马尔可夫链蒙特卡罗(MCMC)的贝叶斯系统发育推断(Yang and Rannala, 1997, Mau et al., 1999)涉及到在给定序列数据(D)的情况下,模拟替代模型参数(φ)和系统发育树的联合后端分布。通过将系统发育模型限制为时间树(见图1),并将系统发育似然与聚合先验耦合,种群历史的参数(θ) Nθ(t),也可以通过抽样后验概率分布(Drummond et al., 2002)同时估计:$$ f_{\theta g \phi}(\theta, g, \phi \mid D)=\frac{1}{\operatorname{Pr}(D)} \operatorname{Pr}(D \mid g, \phi) f_G(g \mid \theta) f_{\Theta}(\theta) f_{\Phi}(\phi) $$ 术语Pr(D∣g, φ)通常被称为系统发育似然,是给定时间树g和替代模型参数的数据的概率。它可以通过修剪算法(Felsenstein, 1981)计算,该算法有效地在树的内部节点上对所有祖先序列状态进行相加。可能性的扩展考虑了站点间的异质性(Yang, 1994)。如果时间树g涉及序列的一个异时样本,那么替换参数φ还包括总体替换率μ,这可以从异时数据中估计,这样就可以在日历尺度上估计种群历史。归一化常数Pr(D)也被称为配分函数或边际似然,它的大小提供了模型支持的度量,尽管它的估计需要先进的MCMC技术(例如,热力学积分或跨维MCMC)。 在确定时间树拓扑的先验密度和聚结/发散时间时,聚结模型将发挥作用。聚结提供了一个概率分布fG(g∣θ),条件是种群规模历史的确定性模型Nθ(t)。它的参数(θ)又可以估计为超参数。给定由边图Eg和聚结次数t = {tn = 0,tn−1,…,t2,t1}组成的n个同时期样本的时间树g = {Eg,t},聚结密度为:$$ f_G(g \mid \theta)=\prod_{i=1}^{n-1}\left[N_\theta\left(t_i\right)^{-1} \exp \left(-\int_{t_i}^{t_{i+1}} \frac{\left(\begin{array}{l} n \\ 2 \end{array}\right)}{N_\theta(t)} d t\right)\right] . $$先验分布fΘ(θ)和fΦ(φ)通常从标准的单变量或多变量分布中选择。
在前一节中,我们已经看到,可以使用系统发育学直接从带有时间标记的基因组数据推断疾病暴发的日期、其源人群和病毒传播历史。系统发育模型主要解决进化史的问题,而动力学模型经常被用来预测未来。预测模型很重要,因为它们提供了预测新流行病结果的某些方面的可能性,并提供了评估大流行病风险和有计划的干预措施的潜在影响的可能性。
系统发育推断是基于遗传数据,如从受感染宿主采样的DNA序列。目前使用这些数据来推断过去信息的模型通常需要简化关于人口规模的假设,例如,保持恒定或服从纯粹的指数增长。另一方面,流行病学家将他们的模型与流行率或发病率数据相匹配。标准的流行病学模型是用一组常微分方程来描述的,这些方程跟踪易感和受感染人数的变化(通常是非线性的)。因此,从生态学的角度来看,系统发育中使用的(受感染个体的)种群大小的简单先验假设是不充分的。
流行病学模型在决定采取何种疾病控制措施以避免或阻止病毒暴发方面发挥着重要作用。通过模型模拟估计隔离、接种疫苗和其他措施的影响,以此作为制定何种公共卫生政策和采取何种行动的决策基础。然而,了解病毒暴发的系统发育历史对于重建传播途径至关重要,这有助于有效管理和未来的预防工作(例如,Cottam等人,2008年)。
决定快速进化RNA病毒多样性的流行病学和生态过程与突变产生的时间尺度相同,并且在种群中是固定的(Holmes, 2004年)。这意味着基因序列数据可以为传播历史提供独立的证据。虽然流行病学数据通常提供关于谁在何时被感染的信息,但它通常不能提供关于传播史的积极证据。因此,这些信息来源的结合应该为更详细的流行病学推断开辟道路,包括接触网络和传播历史的贝叶斯估计(Welch et al., 2011)。
标准的流行病学模型基于宿主间的通量,将宿主群体分为易感(S)、感染(I)和康复或死亡(R)个体。标准模型被称为SI、SIS和SIR。模型的选择基于所考虑疾病的特征、潜伏期的存在、感染后的免疫力等(见专栏2)(Anderson和May, 1991年,Keeling和Rohani, 2008年)。这些模型将焦点限制在每个单元中个体数量的时间演化上,从而掌握了流行病的总体进展。某些疾病特征需要对标准模型进行调整或扩展,例如,纳入无症状感染,在感兴趣的病毒并不总是引起明显症状的情况下,对有症状感染的采样偏向(例如,Aguas等人,2008年)。一个重要的阈值比是基本生殖比R0,即在完全易感人群中由一次原发感染引起的继发感染的预期数量(Diekmann et al., 1990)。流行病学家根据它的价值对疾病的影响进行预测。在经典的确定性流行病学模型中,如果基本生殖比率大于1,则预计会发生流行病。
通常情况下,流行病学家根据经验数据拟合一套合适的确定性微分方程,通常是人群中感染或相关住院人数。因此,该模型可用于估计是否可以通过以下措施来控制流行病:㈠接种疫苗和㈡对易感者进行抗病毒预防、㈢对感染者进行治疗或㈣将感染者与易感者隔离。公共卫生政策的决定往往基于这些估计。
最简单的流行病学模型假设人口中均质混合。在许多情况下,这种假设是不成立的。由于宿主接触动态,病毒感染很容易在学校、城市和农场等社会单位内传播,而在它们之间则不太容易。因此,人口结构一体化是必不可少的。然而,即使在亚种群内,个体动态也可能随机不同(见图4)。这种随机性是可以解释的
近年来,贝叶斯系统发育推断导致了对快速进化病毒的分析爆发。虽然这一爆炸式发展在阐明当前感染因子多样性背后的起源、传播途径和进化速率的多种变化方面卓有成效,但有一个新兴领域有望通过系统发育学和数学流行病学的统一,扩展分子序列数据的概念范围。这一新的系统动力学领域既包括使用系统发育学对经典流行病学参数的推断,也包括旨在研究进化(突变、漂移、达尔文选择)和生态(种群动态和生态随机)过程之间不可避免的相互作用的后果的令人兴奋的新方法。正在进行的研究对进化生物学和分子生态学有更广泛的影响。当一个种群含有不同内在动力特性(如毒性、传播率、恢复率)的基因型时,这种进化和生态的相互作用就会发生。尽管这一条件几乎总是在真实的人群中得到满足,而且往往在形成结果方面起决定性作用,但在流行病学模型中对达尔文选择的数学和理论分析是新兴的系统动力学领域中最具挑战性和研究最少的领域。因此,未来的研究时机已经成熟。
同时,系统动力学研究很可能会迅速发展出统计系统地理学和结构化种群动力学的新方法。
参考文献